include("TKTD_04_boilerplate_dl_ni.jl"); # set up simulation and load default parameters for d magna
Activating project at `c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations`
WARNING: using ProgressMeter.update! in module Main conflicts with an existing identifier.
WARNING: redefinition of constant EMPTVEC. This may fail, cause incorrect answers, or produce other errors.
WARNING: redefinition of constant EMPTVEC. This may fail, cause incorrect answers, or produce other errors.
WARNING: redefinition of constant EMPTVEC. This may fail, cause incorrect answers, or produce other errors.
WARNING: redefinition of constant EMPTVEC. This may fail, cause incorrect answers, or produce other errors.
# using the longispina AMP parameters as default
default_flea_params = copy(default_fleas[:DL_AMP])
# turn off ageing and adaptive plasticity by default
default_flea_params[:h_a_accel] = 0.
default_flea_params[:e_50_adpt] = 1e-10
default_flea_params[:beta_e_adpt] = 10.
# phytoplankton (food) parameters
phyto_params = [copy(default_phyto.Rsub)]
1-element Vector{OrderedDict{Symbol, Any}}:
OrderedDict(:species => "Rsub", :mu_max => 1.78, :q_max => 0.018, :q_min => 0.0011, :q_init => 1.0, :v_max => 0.062, :k_s => 0.0625, :m_max => 0.0, :I_opt => 120.0, :T_opt => 28.0…)
The following code has the purpose to validate the DEB-IBM implementation for each submodel.
An abbreviated version will be provided in the TRACE documentation.
To facilitate comparison of model output with analytical solutions, some changes to the simulation settings are made:
default_global_params[:A_inoculate] = [100.] # => f=1
default_global_params[:renew_medium] = Int64[] # no renewals
default_global_params[:remove_juveniles] = 0:7 |> collect; # remove juveniles immediately
default_global_params[:age_init_max] = 1e-7; # start directly at birth
default_global_params[:data_recording_interval] = 1.;
According To Kooijman (2010), age at birth (i.e., embryonal development time) can be predicted analytically when $g$ is small and $\dot{k}_M$ is large:
$$ a_b = \frac{3l_b}{\dot{k}_M g}\left( 1 + \frac{g}{f} \right) $$, with $L_{max}=\kappa \frac{\left\{ \dot{p}_{A,max} \right\}}{\left[ \dot{p}_M \right]}$, $g=\frac{ [E_G]\dot{v} }{ \kappa \left\{ \dot{p}_{A,max} \right\} }$ and $\dot{k}_M=\frac{\left[ \dot{p}_M \right]}{E_G}$.
This equation also requires that $f$ is equal to scaled reserve density at birth, $e_b$.
$l$ is length $L$ scaled by maximum length $L_{max}$.
We know all the inputs except for $l_b$. So we cannot calculate $a_b$ directly, but we can check for consistency by inserting $l_b$ from the simulation output table, calculating $a_b$ correspondingly, and comparing with $a_b$ in the simulation output table (which is the actual age at birth that results from simulating embryonal development).
function L_max(
kappa::Float64,
p_A_m::Float64,
p_M_V::Float64
)
return kappa*(p_A_m/p_M_V)
end
function g(
E_G::Float64,
v::Float64,
kappa::Float64,
p_A_m::Float64)
return (E_G*v)/(kappa*p_A_m)
end
function k_M(
p_M_V::Float64,
E_G::Float64
)
return p_M_V/E_G
end
k_M (generic function with 1 method)
"""
Calculate age at birth for small g and large k_M.
$(TYPEDSIGNATURES)
"""
function age_at_birth(
L_b::Float64,
L_max::Float64,
k_M::Float64,
g::Float64,
f::Float64
)
let l_b=L_b/L_max
return ((3l_b)/(k_M*g))*(1+(g/f))
end
end
age_at_birth
The code below will run simulation of a life-table experiment for a single test animal, and show output for length and age at birth for a single individual.
flea_params = deepcopy(default_flea_params)
out = life_table(default_global_params, phyto_params, flea_params; n_reps=1)
L_b_sim, a_b_sim = @subset(out.repro, :t_day.==minimum(:t_day))[1,[:L, :age]]
DataFrameRow (2 columns)
| L | age | |
|---|---|---|
| Any | Any | |
| 1 | 0.0172716 | 3.25 |
This code cell will show the analytical solution, making use of the simulated $L_b$. The value that is printed should be equal to "age" in the table above.
L_max_anl = L_max(flea_params[:kappa_0], flea_params[:p_A_m_0], Float64(flea_params[:p_M_V_0]))
g_anl = g(Float64(flea_params[:E_G_0]), flea_params[:v], flea_params[:kappa_0], flea_params[:p_A_m_0])
k_M_anl = k_M(Float64(flea_params[:p_M_V_0]), Float64(flea_params[:E_G_0]))
a_b_anl = age_at_birth(L_b_sim, L_max_anl, k_M_anl, g_anl, 1.)
3.266038320438166
The values are fairly close, but it might be useful to generate repeated predictions and also use different input values for $\kappa$ and $\dot{v}$ to get a better idea.
"""
Get simulated and analytically calculated age at bith, given values for kappa and v.
$(TYPEDSIGNATURES)
"""
function evaluate_a_b(kappa::Float64, v::Float64)
flea_params = deepcopy(default_flea_params)
# smallest E_G_0 value that doesn't break the model (k_M becomes large as a result of this)
flea_params[:kappa_0] = kappa
flea_params[:v] = v
out = life_table(default_global_params, phyto_params, flea_params; n_reps=1)
L_b_sim, a_b_sim = @subset(out.repro, :t_day.==minimum(:t_day))[1,[:L, :age]]
L_max_anl = L_max(flea_params[:kappa_0], flea_params[:p_A_m_0], flea_params[:p_M_V_0])
g_anl = g(flea_params[:E_G_0], flea_params[:v], flea_params[:kappa_0], flea_params[:p_A_m_0])
k_M_anl = k_M(flea_params[:p_M_V_0], flea_params[:E_G_0])
a_b_anl = age_at_birth(L_b_sim, L_max_anl, k_M_anl, g_anl, 1.)
return a_b_sim, a_b_anl
end
vec_a_b_sim = []
vec_a_b_anl = []
kappa_sim = []
for kappa in range(0.1, 0.9, length=10)
try
a_b_sim, a_b_anl = evaluate_a_b(kappa, default_flea_params[:v])
push!(vec_a_b_sim, a_b_sim)
push!(vec_a_b_anl, a_b_anl)
push!(kappa_sim, kappa)
catch
nothing
end
end
p = scatter(
kappa_sim, vec_a_b_anl,
xlabel=L"\kappa\ (-)", ylabel="Embryonal \n development time (d)",
leg=:outertopright, label="Analytical solution",
size=(800,250), leftmargin=5mm, bottommargin=5mm
)
scatter!(kappa_sim, vec_a_b_sim, label="Simulation output")
savefig(plot(p, dpi=300), "figures/implementation_verification_embryonal_development_time.png")
p
vec_a_b_sim = []
vec_a_b_anl = []
v_sim = []
for v in range(0.0022, 0.02, length=10)
try
a_b_sim, a_b_anl = evaluate_a_b(default_flea_params[:kappa_0], v)
push!(vec_a_b_sim, a_b_sim)
push!(vec_a_b_anl, a_b_anl)
push!(v_sim, v)
catch
nothing
end
end
scatter(
v_sim, vec_a_b_anl,
xlabel=L"\dot{v}\ (cm\ d^{-1})", ylabel="Age at birth (d)",
leg=:outertopright, label="Analytical solution",
size=(800,250), leftmargin=5mm, bottommargin=5mm
)
scatter!(v_sim, vec_a_b_sim, label="Simulation output")
The simulation output follows the same pattern as the analytical solution. There is a consistent and very marginal difference in the absolute values of $a_b$. It is not entrirely clear what causes it, but the differences are probably small enough to be irrelevant, and occur for parameter combinations which are biologically irrelevant.
Maximum length (at constant and abundant food) can be directly calculated from parameters: $$ L_{max} = \kappa \frac{ \left\{ \dot{p}_{A,max} \right\} }{ \left[ \dot{p}_M \right] } $$
To verify that $L_{max}$ is replicated correctly, we need to set assimilation efficiency constant and calculate it internally from maximum ingestion and assimilation rates.
flea_params = deepcopy(default_flea_params)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
# increase t_max to be sure that the analytical L_max will be reached
global_params[:t_max] = 56;
# calculate L max analytically
L_max_anl = L_max(flea_params[:kappa_0], flea_params[:p_A_m_0], flea_params[:p_M_V_0])
0.054705986990994986
out = life_table(global_params, phyto_params, flea_params; n_reps=1)
p = @df out.growth plot(
:t_day, :L,
xlabel="Time since birth (d)", ylabel="Structural length (cm)",
leg=:bottomright, label="Simulation output",
ylim=(0, 0.06)
)
hline!([L_max_anl], label="Analytically predicted maximum length")
savefig(plot(p, dpi=300), "figures/implementation_verification_maximum_length.png")
p
flea_params = deepcopy(default_flea_params)
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:e_50_adpt] = 0.75
flea_params[:beta_e_adpt] = 2.0
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.1]
#global_params[:renew_medium] = [0, 2, 5]
1-element Vector{Float64}:
0.1
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
We simulated food being added daily, but never being removed by medium renewals. Since a single animal cannot ingest that much, food accumulates in the medium:
@df out.repro plot(
plot(:t_day, :X_p , linetype=:steppost), size=(400, 200), xlabel="Time (d)", ylabel="Food density (J/L)",
xticks=0:7:21
)
We can use this output to check the shape of the scaled functional response, which should perfectly follow the Type II functional response equation.
This can be compared with the functional response plotted on top:
# calulating the functional response directly
xvals = range(extrema(out.repro.X_p)..., length=100)
f_anl = @. xvals / ((flea_params[:J_X_A_m_0]/flea_params[:F_m_0])+xvals)
# comparison with simulation output
plot(xvals, f_anl, xlabel="Food density (J/L)", ylabel="f", leg=:bottomright, label="Analytical solution", size=(400,200))
@df out.repro scatter!(:X_p, :f, label="Simulation output")
This shows that the functional response is implemented correctly, but it might be useful to check that the scaling by surface area is also correct. We can effectively fix $f$ to a value of 1 by using a very high food density, then plot ingestion rate against surface area and should see a perfectly linear relationship:
flea_params = deepcopy(default_flea_params)
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [1000.]
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
@df out.repro scatter(:L.^2, :J_X, xlabel=L"L^2 (cm^2)", ylabel=L"\dot{J_X} (J/d)", size=(400,300))
Lastly, we need to check that the external food density is updated correctly. We do so by calculating the differential along the vector of output food densities and comparing this with the output ingestion rate, which, in the absence of medium renewals, should be mirror images.
They might not be perfectly reciprocal because neonates are only removed daily and will also ingest some food, but the difference should not exceed a few percent of the total change in food density.
For this, we turn off daily addition of food, so food densities should decline continuously.
flea_params = deepcopy(default_flea_params)
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1/24
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.4]
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
┌ Warning: Inference of survival will currently fail for data_recording_interval !=1. └ @ Main c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations\src\Simulation.jl:75
For adaptive plasticity, there are no analytical solutions that would be useful for implementation verification. We do a quick check that $\kappa$ increases with increasing $e$ and approaches $kappa_0$.
flea_params = deepcopy(default_flea_params)
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:e_50_adpt] = 0.5
flea_params[:beta_e_adpt] = 4.0
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.1]
#global_params[:renew_medium] = [0, 2, 5]
1-element Vector{Float64}:
0.1
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
Note that it is not guaranteed that $\kappa_0$ will be reached. This depends on the value of $\beta_{e,adpt}$.
@df out.repro plot(:e, :kappa, xlim=(0,1), ylim=(0,1), label=L"\kappa_e", xlabel="e", ylabel=L"\kappa", marker=true)
hline!([flea_params[:kappa_0]], label=L"\kappa_0")
Reproduction should occur in discrete intervals, dictated by the time between molts. Note that the plot below indicates reproduction at the time-points when neonates hatch (after embryonal development has been completed), not at the time-points of oviposition.
flea_params = deepcopy(default_fleas.DM_AMP) # use d magna parameters
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.4]
global_params[:data_recording_interval] = 1/24
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
@df out.repro plot(:t_day, :cum_repro, linetype=:steppost, xlabel="Time (d)", ylabel="Cumulative reproduction")
┌ Warning: Inference of survival will currently fail for data_recording_interval !=1. └ @ Main c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations\src\Simulation.jl:75
If we increase $\tau_{molt}$, reproduction should occur less frequently. ($t_{max}$ is increased to demonstrate this).
flea_params = deepcopy(default_fleas.DM_AMP) # use d magna parameters
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:tau_molt] = 7
flea_params[:h_a] = 0.
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.4]
global_params[:data_recording_interval] = 1/24
global_params[:t_max] = 56
out = life_table(global_params, phyto_params, flea_params; n_reps=1);
@df out.repro plot(:t_day, :cum_repro, linetype=:steppost, xlabel="Time (d)", ylabel="Cumulative reproduction")
┌ Warning: Inference of survival will currently fail for data_recording_interval !=1. └ @ Main c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations\src\Simulation.jl:75
The energy in the reproduction buffer should increase continously during adulthood, and decrease in a discrete manner.
flea_params = deepcopy(default_fleas.DM_AMP) # use d magna parameters
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:A_inoculate] = [0.4]
global_params[:data_recording_interval] = 1/24
out = life_table(global_params, phyto_params, flea_params; n_reps=1)
@df out.repro plot(:t_day, :E_R, linetype=:steppost, xlabel="Time (d)", ylabel=L"E_R\ (J)")
┌ Warning: Inference of survival will currently fail for data_recording_interval !=1. └ @ Main c:\Users\shansul\Documents\Research\COM-MET-MIX\WP2\TKTD_calibrations\TKTD_calibrations\src\Simulation.jl:75
If we zoom in to the lowest section of the y-axis, we should see that the reproduction buffer is not emptied entirely, but a small rest remains after oviposition.
This is the amount of energy that was not sufficient to produce yet another egg, and is therefore also not spent.
@df out.repro plot(:t_day, :E_R, linetype=:steppost, xlabel="Time (d)", ylabel=L"E_R\ (J)", ylim=(0,0.1))
The figure below shows the implemented relationship between scaled reserve density and daily survival probability.
x = range(0,1, length=100)
plot(x, exp.(log.(x).*1), xlabel="Scaled reserve density", ylabel="Daily survival probability")
plot!(x, @. exp(log(x)*0.25))
plot!(x, @. exp(log(x)*4))
Survival is simulated at two food levels for three different values of the starvation parameter $\beta_e$.
default_flea_params[:beta_e]
0
flea_params = deepcopy(default_flea_params)
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.4]
global_params[:t_max] = 56
outs = LifeTableDataset[]
e_starve_sweep = [0.125, 0.25, 0.5]
for e_starve in e_starve_sweep
for A_0 in [0.1, 0.4]
global_params[:A_0] = [A_0]
flea_params[:e_starve] = e_starve
flea_params[:beta_e] = 0.5
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.survival[!,:e_starve] .= e_starve
out.survival.food_level .= A_0
push!(outs, out)
end
end
out = concat(outs)
typeof(out)
LifeTableDataset
plot(
(@df @subset(out.survival, :food_level.==0.4) plot(
:t_day, :survival, group=:e_starve, linetype=:steppost,
xlabel=L"Time\ since\ birth\ (d)", ylabel=L"Survival\ rate\ (-)",
legend=false,
leftmargin=5mm, bottommargin=5mm, xticks=0:14:56, title="High food \n (0.4 mg dm/day)"
)),
(@df @subset(out.survival, :food_level.==0.1) plot(
:t_day, :survival, group=:e_starve, linetype=:steppost,
xlabel=L"Time\ since\ birth\ (d)", ylabel=L"Survival\ rate\ (-)",
legend_title="Starvation \n parameter (-)",
leftmargin=5mm, bottommargin=5mm, xticks=0:14:56, title="Low food \n (0.1 mg dm/day)",
titlefontsize=12, legend=:outertopright
)),
lw=1.5, background=:transparent, labelfontsize=14, layout=@layout([a{0.4w} b{0.6w}]),
size=(1200,450), topmargin=5mm, ylim=(0,1)
)
flea_params = deepcopy(default_flea_params)
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:e_starve] = 0.25
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.1]
global_params[:t_max] = 56
outs = LifeTableDataset[]
beta_e_sweep = [0.01, 0.1, 2.0]
for beta_e in beta_e_sweep
flea_params[:beta_e] = beta_e
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.survival.beta_e .= beta_e
push!(outs, out)
end
out = concat(outs)
typeof(out)
LifeTableDataset
@df out.survival plot(
:t_day, :survival, group=:beta_e,
legend=:outertopright, linetype=:steppost,
xlabel=L"Time\ since\ birth\ (d)", ylabel=L"Survival\ rate\ (-)",
legend_title=L"Starvation parameter (-)", size=(800,300),
leftmargin=5mm, bottommargin=5mm, xticks=0:14:56
)
The ageing submodel is checked for plausability by altering the ageing acceleration $\ddot{h}_a$ and comparing survival rates.
flea_params = deepcopy(default_flea_params)
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.4]
global_params[:t_max] = 56
outs = LifeTableDataset[]
h_a_sweep = [0, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1]
for h_a in h_a_sweep
flea_params[:h_a_accel] = h_a
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.survival.h_a .= h_a
push!(outs, out)
end
out = concat(outs)
typeof(out)
LifeTableDataset
As we would expect, high values for $\ddot{h}_a$ cause individuals to die from ageing at a youger age.
@df out.survival plot(
:t_day, :survival, group=:h_a,
legend=:outertopright, linetype=:steppost,
xlabel=L"Time\ since\ birth\ (d)", ylabel=L"Survival\ rate\ (-)",
legend_title=L"Ageing\ acceleration\ (1/d^2)", size=(800,300),
leftmargin=5mm, bottommargin=5mm, xticks=0:14:56
)
TK is controlled by a single parameter $\dot{k}_d$.
The TK module is checked for plausibility by turning off TD (setting all parameters to 0.) and only simulating TK.
Comparison with analytical solutions is possible, but currently considered too time-intesive for the purpose of validating this implementation.
The plot below shows that, for a $k_d$ value of 1 $d^{-1}$, the substance is taken up very fast.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.]
global_params[:A_0] = [0.4]
global_params[:t_max] = 56.
flea_params = deepcopy(default_fleas.DM_AMP)
#flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:TKTD][1].TK[1].params = [1.]
outs = LifeTableDataset[]
C_sweep = [0., 10., 100.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
@df out.repro plot(:t_day, [c[1] for c in :C_int_1], group=:C, linetype=:steppost, leg=:outertopright,
xlabel="Time (d)", ylabel="Scaled damage", legend_title="External concentration",
size=(800,300), bottommargin=5mm, leftmargin=5mm, lw=1.5
)
Lowering $k_d$ to 0.1 leads, expectedly, to a more gradual uptake.
flea_params[:TKTD][1].TK[1].params = [.1]
outs = LifeTableDataset[]
C_sweep = [0., 10., 100.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
@df out.repro plot(:t_day, [c[1] for c in :C_int_1], group=:C, linetype=:steppost, leg=:outertopright,
xlabel="Time (d)", ylabel="Scaled damage", legend_title="External concentration",
size=(800,300), bottommargin=5mm, leftmargin=5mm, lw=1.5
)
flea_params[:TKTD][1].TK[1].params = [.05]
outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
@df out.repro plot(:t_day, [c[1] for c in :C_int_1], group=:C, leg=:outertopright,
xlabel="Time (d)", ylabel="Scaled damage", legend_title="External concentration",
size=(800,300), bottommargin=5mm, leftmargin=5mm, lw=1.5, linetype=:steppost
)
outs = LifeTableDataset[]
kd_sweep = 10 .^ range(log10(0.005), log10.(1), length=100)
for k_d in kd_sweep
global_params[:C][1] = 100.
flea_params[:TKTD][1].TK[1].params = [k_d]
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.k_d .= k_d
push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
@df @subset(out.repro, round.(:t_day, sigdigits=2).==21) plot(:k_d, [c[1] for c in :C_int_1], xscale=:log10, ylim=(0,101),
xlabel=L"\dot{k}_d\ (d^{-1})", ylabel="Concentration", label="Scaled damage at day 21"
)
hline!([100.], leg=:bottomright, label="External concentration")
TK models can be linked for different PMoAs, or independent of each other.
When link_TK=true, the damage values for all PMoAs should be identical.
flea_params[:TKTD][1].TK[1].params = [.05]
flea_params[:TKTD][1].link_TK = true
outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day)
# matrix of damage values per TK component
damage_matrix = vcat(out.repro.C_int_1...)
# this statement needs to be true
damage_matrix[:,1]==damage_matrix[:,2]==damage_matrix[:,3]==damage_matrix[:,4]==damage_matrix[:,5]
true
When link=false, damage dynamics are different for every PMoA.
flea_params[:TKTD][CU].TK[1].params = [.05]
flea_params[:TKTD][CU].TK[2].params = [1.0]
flea_params[:TKTD][CU].TK[3].params = [0.]
flea_params[:TKTD][CU].link_TK = false
outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.repro, :t_day);
The figure below should show slow TK for TK component 1, fast TK for TK component 2, and no increase in damage for TK component 3.
@df out.repro plot(
groupedlineplot(:t_day, [c[1] for c in :C_int_1], :C, marker=true, title="TK component 1", ylim=(0,100)),
groupedlineplot(:t_day, [c[2] for c in :C_int_1], :C, marker=true, title="TK component 2"),
groupedlineplot(:t_day, [c[3] for c in :C_int_1], :C, marker=true, title="TK component 3", ylim=(0,100)),
xlabel="Time (d)", ylabel="Scaled \n damage", bottommargin=5mm, leftmargin=5mm, size=(1000,600),
background_color=:transparent
)
For the TD module, we first check for the presence/absence of effects on reproduction, growth and survival for each mode of action. Reproduction should be affect by all modes of action. Growth should be affected by all by reproduction efficiency. Maintenance costs and Assimilation efficiency should affect survival by triggering starvation mortality
Top row shows carapace length (cm) over time, bottom row cumulative reproduction. Each column is an exposure concentration.
outs = LifeTableDataset[]
for C in C_sweep
global_params[:C][CU] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro[!,:C] .= C
out.growth[!,:C] .= C
out.survival[!,:C] .= C
push!(outs, out)
end
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.
flea_params[:TKTD][CU].TD = default_flea_params[:TKTD][1].TD[[G]]
flea_params[:TKTD][CU].TD[1].DRC.params = (100., 2., 5.)
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true
outs = LifeTableDataset[]
C_sweep = [0.01, 10., 100., 1000.]
for C in C_sweep
global_params[:C][CU] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro[!,:C] .= C
out.growth[!,:C] .= C
out.survival[!,:C] .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))
combine(groupby(out.growth, :C)) do df
L0 = df.carapace_length[1]
L7 = df[argmin(abs.(df.t_day .- 7)),:carapace_length]
return DataFrame(
SGR = (log(L7)-log(L0))/7
)
end |> x-> plot(x.C, x.SGR, xlabel="External concentration", ylabel="7-day somatic growth rate", xscale=:log10, lw=2, title="Effects on somatic growth rate", marker=true)
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.0
flea_params[:TKTD][CU].TD = default_flea_params[:TKTD][1].TD[[M]]
flea_params[:TKTD][CU].TD[1].DRC.params = (100., 2., 5.)
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true
outs = LifeTableDataset[]
C_sweep = [0.01, 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))
combine(groupby(out.growth, :C)) do df
L0 = df.carapace_length[1]
L7 = df[argmin(abs.(df.t_day .- 7)),:carapace_length]
return DataFrame(
SGR = (log(L7)-log(L0))/7
)
end |> x-> plot(x.C, x.SGR, xlabel="External concentration", ylabel="7-day somatic growth rate", xscale=:log10, lw=2, title="Effects on somatic growth rate", marker=true)
Assimilation efficiency causes effects on all variables, including survival if $\dot{h}_{s,e}$ is large enough.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[A].DRC.params = (100., 2., 1.)
outs = LifeTableDataset[]
C_sweep = [0.01, 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[R].DRC.params = (100., 2., 1.)
outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))
Same for slow TK to check that the connection between TK and TD works as desired.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [.01]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[EO].DRC = DRCModel(
LL3GM,
(100., 2., 1.)
)
outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))
Increase in $\kappa$ should result in increased growth including maximum length and decreased reproduction.
default(leg=false)
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:e_50_adpt] = 1e-3
flea_params[:beta_e_adpt] = 10.
flea_params[:TKTD][CU].TK[KP].params = [1.]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[KP].DRC = DRCModel(
LL3GM,
(100., 2., 1.)
)
outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450), bottommargin=5mm)
Decrease in $\kappa$ should result in decreased growth and (slightly) increased reproduction.
The stimulating effect on reproduction should be relatively subtle compared to the stimulating effect on growth resulting from increased $\kappa$,
because it implies a gross decrease in growth, therefore a gross decrease in reproduction, as well as a gross increase in reproduction.
The parameter combination below should lead to a net decrease in reproduction.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [1.]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[KP].DRC = DRCModel(
LL3GM,
(100., 2., -.5)
)
outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))
In the code below, the TK parameter is lowered which should produce a hormetic effect on reproduction.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [.1]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[KP].DRC = DRCModel(
LL3GM,
(100., 2., -.5)
)
outs = LifeTableDataset[]
C_sweep = [0., 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,120))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450))
outs = LifeTableDataset[]
C_sweep = geomrange(0.1, 1000.)
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
LifeTableDataset(11000×60 DataFrame Row │ t_day unique_id species_idx kappa kappa_EX q_p_rel J_X ⋯ │ Any Any Any Any Any Any Any ⋯ ───────┼──────────────────────────────────────────────────────────────────────── 1 │ 0.0 1 1 0.66472 1.0 1.0 0.222054 ⋯ 2 │ 1.0 1 1 0.672075 1.0 1.0 0.459321 3 │ 2.0 1 1 0.668482 1.0 1.0 0.758668 4 │ 3.0 1 1 0.667277 1.0 1.0 1.10867 5 │ 4.0 1 1 0.666679 1.0 1.0 1.49666 ⋯ 6 │ 5.0 1 1 0.666323 1.0 1.0 1.91205 7 │ 6.0 1 1 0.666089 1.0 1.0 2.34598 8 │ 7.0 1 1 0.665924 1.0 1.0 2.79109 ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ 10994 │ 15.0 1 1 0.33804 1.0 1.0 2.16702 ⋯ 10995 │ 16.0 1 1 0.337747 1.0 1.0 2.24583 10996 │ 17.0 1 1 0.337501 1.0 1.0 2.31961 10997 │ 18.0 1 1 0.337291 1.0 1.0 2.3886 10998 │ 19.0 1 1 0.337111 1.0 1.0 2.45305 ⋯ 10999 │ 20.0 1 1 0.336957 1.0 1.0 2.51317 11000 │ 21.0 1 1 0.336823 1.0 1.0 2.56923 53 columns and 10985 rows omitted, 11000×5 DataFrame Row │ t_day rep L carapace_length C │ Any Int64 Any Float64 Float64 ───────┼─────────────────────────────────────────────────── 1 │ 0.0 1 0.0264072 0.11949 0.1 2 │ 1.0 1 0.037256 0.168579 0.1 3 │ 2.0 1 0.047536 0.215095 0.1 4 │ 3.0 1 0.0572316 0.258967 0.1 5 │ 4.0 1 0.0663241 0.300109 0.1 6 │ 5.0 1 0.0748312 0.338603 0.1 7 │ 6.0 1 0.0827806 0.374573 0.1 8 │ 7.0 1 0.0902034 0.40816 0.1 ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ 10994 │ 15.0 10 0.0791303 0.358056 1000.0 10995 │ 16.0 10 0.0805421 0.364444 1000.0 10996 │ 17.0 10 0.0818418 0.370325 1000.0 10997 │ 18.0 10 0.0830385 0.37574 1000.0 10998 │ 19.0 10 0.0841408 0.380727 1000.0 10999 │ 20.0 10 0.0851562 0.385322 1000.0 11000 │ 21.0 10 0.0860918 0.389556 1000.0 10985 rows omitted, 1100×4 DataFrame Row │ t_day num_surviving survival C │ Float64 Int64 Float64 Float64 ──────┼─────────────────────────────────────────── 1 │ 0.0 10 1.0 0.1 2 │ 1.0 10 1.0 0.1 3 │ 2.0 10 1.0 0.1 4 │ 3.0 10 1.0 0.1 5 │ 4.0 10 1.0 0.1 6 │ 5.0 10 1.0 0.1 7 │ 6.0 10 1.0 0.1 8 │ 7.0 10 1.0 0.1 ⋮ │ ⋮ ⋮ ⋮ ⋮ 1094 │ 15.0 10 1.0 1000.0 1095 │ 16.0 10 1.0 1000.0 1096 │ 17.0 10 1.0 1000.0 1097 │ 18.0 10 1.0 1000.0 1098 │ 19.0 10 1.0 1000.0 1099 │ 20.0 10 1.0 1000.0 1100 │ 21.0 10 1.0 1000.0 1085 rows omitted)
@df @subset(out.repro, :t_day.==maximum(:t_day)) lineplot(
:C, :cum_repro, xscale=:log10, marker=true, xlabel="External concentration", ylabel="21-day reproduction",
title=L"PMoA $KP$: Effect on 21-day reproduction"
)
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [.1]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[EO].DRC = DRCModel(
LL3AR,
(100.0, 2.0, 0.5)
)
outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(:t_day, :cum_repro, group=:C, layout=(1,4), ylim=(0,200))
ps = @df out.survival plot(:t_day, :survival, group=:C, layout=(1,4), ylim=(0,1.01))
plot(pg, pr, ps, layout=(3,1), xticks=0:7:21, xlabel="Time (d)", lw=2, size=(1000, 450), bottommargin=5mm)
outs = LifeTableDataset[]
C_sweep = geomrange(0.1, 1000.)
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
LifeTableDataset(11000×60 DataFrame Row │ t_day unique_id species_idx kappa kappa_EX q_p_rel J_X ⋯ │ Any Any Any Any Any Any Any ⋯ ───────┼──────────────────────────────────────────────────────────────────────── 1 │ 0.0 1 1 0.66472 1.0 1.0 0.222054 ⋯ 2 │ 1.0 1 1 0.672075 1.0 1.0 0.459321 3 │ 2.0 1 1 0.668482 1.0 1.0 0.758668 4 │ 3.0 1 1 0.667277 1.0 1.0 1.10867 5 │ 4.0 1 1 0.666679 1.0 1.0 1.49666 ⋯ 6 │ 5.0 1 1 0.666323 1.0 1.0 1.91205 7 │ 6.0 1 1 0.666089 1.0 1.0 2.34598 8 │ 7.0 1 1 0.665924 1.0 1.0 2.79109 ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱ 10994 │ 15.0 1 1 0.665441 1.0 1.0 6.21733 ⋯ 10995 │ 16.0 1 1 0.665421 1.0 1.0 6.59176 10996 │ 17.0 1 1 0.665405 1.0 1.0 6.95057 10997 │ 18.0 1 1 0.665391 1.0 1.0 7.29334 10998 │ 19.0 1 1 0.665382 1.0 1.0 7.62018 ⋯ 10999 │ 20.0 1 1 0.665374 1.0 1.0 7.9311 11000 │ 21.0 1 1 0.665369 1.0 1.0 8.22618 53 columns and 10985 rows omitted, 11000×5 DataFrame Row │ t_day rep L carapace_length C │ Any Int64 Any Float64 Float64 ───────┼─────────────────────────────────────────────────── 1 │ 0.0 1 0.0264072 0.11949 0.1 2 │ 1.0 1 0.037256 0.168579 0.1 3 │ 2.0 1 0.047536 0.215095 0.1 4 │ 3.0 1 0.0572316 0.258967 0.1 5 │ 4.0 1 0.0663241 0.300109 0.1 6 │ 5.0 1 0.0748312 0.338603 0.1 7 │ 6.0 1 0.0827806 0.374573 0.1 8 │ 7.0 1 0.0902034 0.40816 0.1 ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ 10994 │ 15.0 10 0.134172 0.607113 1000.0 10995 │ 16.0 10 0.138129 0.625018 1000.0 10996 │ 17.0 10 0.141816 0.641704 1000.0 10997 │ 18.0 10 0.145253 0.657253 1000.0 10998 │ 19.0 10 0.148455 0.671742 1000.0 10999 │ 20.0 10 0.151439 0.685242 1000.0 11000 │ 21.0 10 0.154218 0.697821 1000.0 10985 rows omitted, 1100×4 DataFrame Row │ t_day num_surviving survival C │ Float64 Int64 Float64 Float64 ──────┼─────────────────────────────────────────── 1 │ 0.0 10 1.0 0.1 2 │ 1.0 10 1.0 0.1 3 │ 2.0 10 1.0 0.1 4 │ 3.0 10 1.0 0.1 5 │ 4.0 10 1.0 0.1 6 │ 5.0 10 1.0 0.1 7 │ 6.0 10 1.0 0.1 8 │ 7.0 10 1.0 0.1 ⋮ │ ⋮ ⋮ ⋮ ⋮ 1094 │ 15.0 10 1.0 1000.0 1095 │ 16.0 10 1.0 1000.0 1096 │ 17.0 10 1.0 1000.0 1097 │ 18.0 10 1.0 1000.0 1098 │ 19.0 10 1.0 1000.0 1099 │ 20.0 10 1.0 1000.0 1100 │ 21.0 10 1.0 1000.0 1085 rows omitted)
PMoA Embryonal investement leads to an increase in cumulative reproduction.
p = @df @subset(out.repro, :t_day.==maximum(:t_day)) plot(
lineplot(
:C, :E_o, xscale=:log10, xlabel="External concentration", ylabel="Energy investment per egg (J)", lw=2
),
lineplot(
:C, :cum_repro ./ 100., xscale=:log10, xlabel="External concentration", ylabel="Control-normalized \n 21-day reproduction",
lw=2,
)
)
savefig(plot(p, dpi=300), "figures/pmoa_E_repro.png")
p
Each embryo starts with less energy, so there should be implications for neonate fitness.
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [100.]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[1].params = [.1]
flea_params[:TKTD][CU].link_TK = true
flea_params[:TKTD][CU].TD[EO].DRC = DRCModel(
LL3AR,
(100.0, 2.0, 0.5)
)
DRCModel(LL3AR, (100.0, 2.0, 0.5))
global_params[:remove_juveniles] = []
Any[]
neonate_fitness = DataFrame()
C_sweep = geomrange(0.1, 1000.)
for C in C_sweep
global_params[:C][1] = C
m = initialize_model(global_params, phyto_params, [flea_params])
run!(m)
all_daphnids = daphnid_output(m, Inf)
all_daphnids[:,:C] .= C
append!(neonate_fitness, all_daphnids)
end
combine(groupby(@subset(neonate_fitness, :unique_id.>1), [:unique_id, :C])) do df
birth = @subset(df, :t_day.==minimum(:t_day))
DataFrame(
e_b = birth.e,
E_b = birth.E,
L_b = birth.L,
a_b = birth.age
)
end |>
x-> @df x plot(
lineplot(:C, :E_b, ylim=(0, maximum(:E_b)), xscale=:log10),
lineplot(:C, :e_b, ylim=(0, maximum(:e_b)), xscale=:log10),
lineplot(:C, :L_b, ylim=(0, maximum(:L_b)), xscale=:log10),
lineplot(:C, :a_b, ylim=(0, maximum(:a_b)), xscale=:log10),
xlabel="Exposure concentration", ylabel=""
)
include("src/DEB.jl")
WARNING: redefinition of constant EMPTVEC. This may fail, cause incorrect answers, or produce other errors.
apply_effects!
flea_params[:TKTD][CU].TD
7-element Vector{TDModel}:
TDModel(DRCModel(LL3GM, (0.0, 0.0, 0.0)), :s_G_z, prod)
TDModel(DRCModel(LL3GM, (0.0, 0.0, 0.0)), :s_M_z, prod)
TDModel(DRCModel(LL3AR, (0.0, 0.0, 0.0)), :s_A_z, prod)
TDModel(DRCModel(LL3AR, (0.0, 0.0, 0.0)), :s_R_z, prod)
TDModel(DRCModel(LL3GUTS, (100.0, 2.0, 1.0)), :h_z, sum)
TDModel(DRCModel(LL3AR, (0.0, 0.0, 0.0)), :s_E_o_z, prod)
TDModel(DRCModel(LL3AR, (0.0, 0.0, 0.0)), :s_kappa_z, prod)
global_params = deepcopy(default_global_params)
global_params[:data_recording_interval] = 1.
global_params[:A_inoculate] = [0.4]
flea_params = deepcopy(default_fleas.DM_AMP)
flea_params[:kappa_EX_min] = flea_params[:kappa_EX_max] = 1. #flea_params[:p_A_m_0] / flea_params[:J_X_A_m_0]
flea_params[:cv] = 0.
flea_params[:h_a_accel] = 0.
flea_params[:beta_e] = 0.1
flea_params[:TKTD][CU].TK[H].params = [1.]
flea_params[:TKTD][CU].link_TK = false
flea_params[:TKTD][CU].TD[H].DRC = DRCModel(
LL3GUTS,
(100.0, 2.0, 1.0)
)
outs = LifeTableDataset[]
C_sweep = [0.1, 10., 100., 1000.]
for C in C_sweep
global_params[:C][1] = C
out = life_table(global_params, phyto_params, flea_params; n_reps=10)
out.repro.C .= C
out.growth.C .= C
out.survival.C .= C
push!(outs, out)
end
out = concat(outs)
sort!(out.growth, :t_day)
sort!(out.repro, :t_day)
pg = @df out.growth plot(
:t_day, :carapace_length, group=:C,
xlim=(-1,22), ylim=(0,0.8), layout=(1,4), title=hcat(C_sweep...),
size=(800,300), lw=1.5
)
pr = @df out.repro plot(
:t_day, :cum_repro, group=:C,
layout=(1,4), xlim=(-1,22), ylim=(0,150)
)
ps = @df out.survival plot(
:t_day, :survival, group=:C,
layout=(1,4), ylim=(0,1.01)
)
plot(
pg, pr, ps,
layout=(3,1), xticks=0:7:21,
xlabel="Time (d)", lw=2, xlim=(-1,22),
size=(1000, 450), bottommargin=5mm
)
An algal sinking rate is implemented, as this can contribute substantially to food density in the environment.
It is checked whether the algal sinking rate leads to expected decrease in algal biomass in the medium/
global_params = copy(default_global.zooplankton.population_renewal)
phyto_params = [copy(default_phyto.Rsub)]
1-element Vector{OrderedDict{Symbol, Any}}:
OrderedDict(:species => "Rsub", :mu_max => 1.78, :q_max => 0.018, :q_min => 0.0011, :q_init => 1.0, :v_max => 0.062, :k_s => 0.0625, :m_max => 0.0, :I_opt => 120.0, :T_opt => 28.0…)